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Interacting  fields  approach  for  evolving  spatial 
phenomena:  application  to  erosion  simulation 
for  optimized  land  use 

We  propose  a  computational  framework  and  strategies  for  performing  tasks  necessary  for  evaluation  and 
optimization  of  land  use  management  within  an  advanced  GIS  modeling  environment.  Such  tasks 
involve  modeling  of  landscape  processes,  simulation  of  impact  of  human  activities  on  these  processes 
and  optimization  of  preventive  measures  aimed  at  creating  sustainable  landscapes.  A  typical  example  of 
interaction  between  landscape  processes  and  human  activity  is  water  erosion,  a  natural  process  which 
can  be  accelerated  or  greatly  reduced  by  human  intervention.  In  our  implementation,  the  underlying  2D 
water  and  sediment  transport  equations  are  solved  within  approximate  diffusive  wave  hydrologic  and 
detachment/transport  capacity  models  using  stochastic  methods  (Monte  Carlo).  The  methods  employ 
distributed  input  parameters  and  enable  us  to  simulate  the  impact  of  complex  terrain,  various  soils  and 
land  cover  changes  on  the  spatial  distribution  of  erosion  and  deposition.  The  presented  approach 
together  with  land  use  optimization  scenario  is  illustrated  by  an  application  to  an  experimental  farm  in 
Germany,  carried  out  within  3D  dynamic  GIS  environment. 
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INTRODUCTION 

Efforts  to  balance  the  economic  development  with  environmental  protection  increase  the  demand  for 
simulation  tools  enabling  predictions  of  the  human  impact  on  landscape.  In  order  to  prevent  irreversible 
changes  and  avoid  costly,  ineffective  solutions  the  simulation  tools  should  provide  detailed  spatial  and 
temporal  distributions  of  modeled  phenomena.  Statistical  averages  for  the  entire  study  areas  or 
predictions  only  for  a  certain  point,  such  as  watershed  outlet  are  often  insufficient.  New  developments  in 
GIS,  especially  support  for  multivariate  temporal  data  processing,  analysis,  and  visualization  (Mitasova 
et  al.  1995,  GRASS4.2  libraries)  make  such  simulations  possible  and  GIS  plays  an  important  role  in  the 
development  and  applications  of  distributed  landscape  process  models  (e.g.,  Engel  1995(review  with 
links),  Vieux  1995,  Saghafian  1995,  Leavesley  et  al.  1995  ). 

The  goal  of  our  paper  is  to  outline  a  methodological  framework  for  distributed  models  based  on  the 
solution  of  the  "first  principles"  master  equations  for  multi-variate  fields  and  use  these  tools  for  the 
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distributed  landuse  scenarios  optimization.  The  basic  platforms  of  our  approach  are  described  as 
follows: 

•  Representation  of  phenomena  as  multivariate  fields  (i.e.  as  genuine  distributed  objects).  This 
advances  lumped  or  semi-lumped  description  into  the  high  resolution  distributed  one  with 
advantages  for  spatial  and  temporal  analysis.  Such  approach  requires  multivariate  tools  which  deal 
efficiently  with  transformation  from  one  discrete  format  to  another,  and  with  the  possibilities  of 
continuous  representation  for  calculation  of  gradients,  curvatures  and  other  quantities.  For  these 
purposes  we  use  the  multivariate  splines  developed  previously  (Mitasova  et  al.  1995)  and  tested  on 
a  variety  of  2D,  3D  and  4D  data  (Table  1.). 

•  Phenomena  description  and  prediction  based  on  solving  master  equations  which  determine 
the  configuration  or  evolution  of  corresponding  fields.  This  enables  us  to  perform  simulations 
which  are  formulated  in  terms  of  fundamental  physical  processes  such  as  flux,  diffusion,  etc.,  and 
we  can  use  the  known  mathematical  and  physical  apparatus  to  analyze  the  solutions.  This  level  of 
rigour  avoids  vague  concepts  which  often  characterize  various  semi-empirical  schemes  and  allows 
us  to  clearly  specify  inputs,  governing  parameters  and  to  understand  the  character  of  processes 
involved  as  well  as  final  outputs.  We  can  also  build  upon  experience  from  other  disciplines  in 
using  efficient  and  robust  methods  for  solving  the  underlying  master  equations. 

•  Formulation  of  cost  functionals  and  their  optimizations.  The  proper  quantification  of  objectives 
is  very  important  for  an  efficient  solving  of  desired  landuse  practices.  It  is  now  well  understood 
that  the  formulation  of  an  appropriate  cost  or  objective  functional  provides  a  powerful  strategy 
how  to  deal  with  this  task.  The  cost  functional  depends  on  the  fields  and  includes  also  various 
conditions  or  restrictions.  The  human  actions  can  change  the  character  or  properties  of  some  of  the 
input  fields  or  modify  their  future  evolution.  The  "space  of  human  actions"  is  then  explored  for 
estimating  the  optimal  solutions. 

We  illustrate  these  three  fundamental  concepts  on  a  case  of  erosion  prevention  by  optimized  landuse 
scenario.  Water  erosion,  with  its  economical  and  ecological  impact,  is  a  typical  example  of  problems 
targeted  by  our  approach.  It  is  a  genuine  space-time  distributed  phenomenon  with  several  natural 
components  such  as  terrain,  soils,  cover  and  climate  effects.  These  can  be  naturally  and  elegantly 
represented  by  multivariate  fields  (terrain  surface,  terrain  cover,  water  and  sediment  distributions,  water 
and  sediment  fluxes,  soil  distribution,  etc).  The  erosion  processes  can  be  described  by  master  continuity 
and  momentum  conservation  equations.  Some  of  these  fields  such  as  land  cover  can  be  influenced  or 
changed  by  human  intervention  and  therefore  will  naturally  enter  the  optimization  process. 

In  the  following  sections  we  describe  the  working  examples  of  advanced  methods  and  approaches 
which  enabled  us  to  start  from  basic  input  data  (terrain,  cover,  soils,  etc)  and  get  to  the  actual  creative 
process  of  landuse  optimization  in  a  fully  distributed  manner. 


LANDSCAPE  CHARACTERIZATION  IN  3D  SPACE  AND  TIME 

The  basic  objects  which  enter  into  distributed  models  of  landscape  processes  are  given  by  functions 
which  depend  on  the  position  in  3D  space  and  time:  multivariate  scalar  and  vector  fields.  These  fields 
represent  various  phenomena  such  as  terrain,  soil  properties,  land  cover,  fluxes  of  matter.  They  are 
usually  represented  in  a  GIS  database  in  a  discrete  form  as  sets  of  points  (sites,  lines  or  polygons)  or 
rasters.  However,  they  can  be  transformed  to  continuous  representations  using  expansions  in  an 


appropriate  basis  set  such  as  multivariate  regularized  spline  with  tension  (Mitasova  et  al.  1995),  as 
illustrated  by  examples  in  Table  1.  or  by  Hargrove  et  al.  (1995).  Such  representation  is  not  restricted  to 
continuous  fields,  an  example  of  effective  handling  of  surfaces  with  faults  using  splines  and  GIS  tools 
was  developed  by  Cox  et  al.  1994.  Phenomena  represented  by  classes,  such  as  vegetation  or  land  use, 
can  be  represented  as  fields  with  faults  as  well  (Table  1.,  land  cover). 

Table  1.  Examples  of  landscape  phenomena  representations  by  multivariate  fields 
(  Click  on  the  image  to  retrieve  full  size  picture  or  animation) 


Representations  of  natural  phenomena  as  multivariate  fields  modeled  by  bi-,  tri-,  and 

quad- variate  splines 


phenomenon  (field) 

3D  (dynamic)  "map" 

GIS  format  (discrete  representation) 

elevation 

points  (x,y,z),  lines  (contours),  2D  raster 
(DEM) 

elevation  gradient  and  curvatures 

2D  rasters  derived  from  z=f(x,y) 

precipitation 

points(x,y,z,p,t),  2D  raster  time  series 

soil  horizons 

points  (x,y,z,w),  2D  raster  vertical  series 

land  cover 

polygon,  raster 

underground  concentrations  of 
chemicals 

points  (x,y,z,t,w),  3D  raster  time  series 

concentration  of  chemicals  in  water 

points  (x,y,z,t,w),  3D  raster  time  series 

The  fields  are  static  or  evolving  in  time.  The  evolution  can  be  monitored  and  data  from  monitoring  can 
be  stored  in  a  GIS  database  in  various  formats,  for  example,  as  time  series  of  sites  which  can  be  further 
transformed  to  time  series  of  2D  or  3D  rasters  (e.g.,  Table  l.,chem.  concentrations).  To  predict  the 
future  states  of  these  fields,  we  need  to  understand  the  processes  controlling  their  evolution  and 
formulate  models  simulating  their  fundamental  behavior.  The  GIS  software  is  being  enhanced  to  support 


such  simulations  by  providing  the  adequate  data  structures,  including  2D  and  3D  floating  point  raster 
data  (Waupotitsch  and  Shapiro  1995),  multidimensional  site  data  (McCauley  1995),  support  for 
temporal  data  (Brown  and  Shapiro  1995),  methods  for  transformation  between  discrete  and  continuous 
representations  using  multivariate  interpolation  by  radial  basis  functions  (Mitasova  et  al.  1995),  and 
tools  for  interactive  multidimensional  dynamic  cartography  (Brown  et  al.  1995). 

FLUXES  IN  LANDSCAPE 

The  problem  of  landscape  process  simulations  related  to  land  use  change  have  been  attacked  by  a  variety 
of  approaches  such  as  rule-based  models,  cellular  automata,  probability  transitions  cell  based-models 
(Berry  et  al.  1995)  or  as  continuous  fields  and  differential  equations  (Maidment  1995).  In  our  approach 
to  modeling  of  erosion  processes  we  start  from  the  continuous  formulation  and  we  use  the  fundamental 
conservation  laws  (matter  and  approximate  momentum  conservation). 

Water  flux 

One  of  the  primary  processes  in  landscape,  influencing  the  distribution  of  soils,  plants  and  people  is 
water  flow.  While  there  are  numerous  empirical  and  process  based  models  for  modeling  the  waterflux  in 
streams  and  rivers,  spatial  distribution  of  water  on  complex  hillslopes  (crucial  for  soils  and  organisms)  is 
still  being  modeled  by  rather  rough  estimates  unable  to  provide  sufficient  detail  for  modeling  of  some 
important  water  related  processes,  such  as  erosion. 

In  general,  the  overland  flow  is  described  by  the  continuity  equation: 

(1) 

q(M)-MMMM)  (2) 

where  h  is  the  water  depth  [m],  i  is  the  rainfall  excess=rainfall-infiltration  [m/s],  q  is  the  water  flux 
[m.m/s],  v  is  the  flow  velocity  [m/s],  r=(x,y)  is  the  position  [m],  and  t  is  the  time.  The  velocity  v  is 
related  to  h  by  Mannings  or  Chezy  law  and  the  momentum  conservation  equation  (Haan  1994). 

Current  GISs  provide  the  tools  for  the  simplest  approximate  solution  of  equation  (1)  for  steady  state 
flow  and  constant  velocity,  based  on  the  upslope  contributing  area.  There  are  numerous  algorithms 
available  for  its  estimation,  such  as  D8  (Figure  la),  or  an  improved  approach  based  on  the  vector-grid 
algorithm  (Mitasova  et  al.  1995,  Figure  lb).  While  these  geometric  approaches  provide  enough 
information  for  a  wide  range  of  applications  (Moore  et  al.  1992),  they  become  problematic  when  applied 
to  areas  with  spatially  variable  cover  and  complex  terrains  with  significant  spatial  variations  in  flow 
velocity. 

A  more  realistic  approximation  which  takes  into  account  the  flow  velocity  is  based  on  the  steady  state 
solution  dr(T)  in  the  2D  kinematic  wave  (Figure  lc)  and  2D  approximate  diffusive  wave  (Figure  Id) 
(Mitas  and  Mitasova,  in  prep.).  The  2D  approximate  diffusive  wave  solution  can  be  found  by  solving  the 
steady  state  of  ( 1)  which  can  be  written  in  the  operator  form  as: 

PF[^3(r5)  -  -|V>[^(r)]  +  V  •  (fc(r)v(r)]  -  *(r)  (3) 

A 

where  e  is  the  diffusion  constant  and  W  is  the  operator.  We  use  a  stochastic  method  to  solve  the  equation 
(3).  This  approach  is  based  on  the  representation  of  the  solution  h  by  a  large  set  of  random  walkers 


(sampling  points)  which  are  propagated  according  to  the  Green’s  function  corresponding  to  the  inverse 
operator  W.  The  diffusion  term  in  (3)  describes  (approximately)  a  backwater  effect  and  also  helps  to 
reduce  the  artificial  features  in  water  surfaces  on  hillslopes  (Figure  la,b,c)  caused  by  flowtracing  on  a 
regular  discrete  grid,  which  make  the  application  of  such  water  surfaces  in  erosion/deposition  models 
problematic.  This  approach  can  be  extended  also  for  non-stationary  event  based  modeling. 


m 

1 0.5000 
10.1000 
1 0.0100 
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Figure  1.  Steady  state  water  depth  estimated  by  a)  upslope  contributing  area  using  D8  algorithm,  b) 
upslope  contributing  area  using  vector-grid  algorithm  c)  2D  kinematic  wave  approximation,  d)  2D 
approximate  diffusive  wave 

Other  approaches  which  solve  for  temporal  and  spatial  distribution  of  water  depth  during  storm  events 
are,  for  example,  kinematic  (Garrote  and  Brass  1995,  Vieux  1995)  or  diffusive  wave  models  (Saghafian 
1995,  Figure  le),  already  integrated  with  GIS  (r.water.fea,  r.hydro.CASC2d).  The  choice  of  hydrology 
model  complexity  and  realism  depends  on  the  type  of  application  and  for  our  current  efforts  steady  state 
provides  adequate  information  for  assessing  the  impact  of  water  flow  on  landscape  at  the  time  scale  of 
days  to  years.  The  steady  state  solutions  are  also  consistent  with  the  new  generation  erosion  model 
Water  Erosion  Prediction  Project  (WEPP)  (Flanagan  and  Nearing  1995). 

Sediment  flux  and  erosion/deposition  in  complex  terrain 

Soil  erosion  involves  detachment,  transport  and  deposition.  The  interaction  between  soil  detachment  and 
sediment  transport  is  controlled  by  water  flux,  terrain,  soil  and  cover.  This  interaction  is  very  difficult  to 
capture  by  traditional  empirical  models  or  models  based  on  the  geometrical  analysis  of  terrain.  While 
some  of  these  models  provide  adequate  tools  for  a  qualitative  assessment  of  erosion  risk  for  large  areas 
with  complex  terrain  (Moore  and  Wilson  1992,  Mitasovaet  al.  1996),  they  are  insufficient  for  modeling 
of  impact  of  spatially  variable  landuse  and  simulation  of  erosion  protection  measures  effectiveness. 

The  basic  relationship  for  fundamental  erosion  processes  is  continuity  of  mass.  For  erosion  by  2D 
overland  flow,  the  continuity  equation  is  (Foster  and  Meyer  1972,  Govindaraju  1991) 

+  y  .  -  sources  -  sinks  (4) 

at 

where  q^s=r_5  c(r,t) q  is  the  sediment  flux  [kg/(ms)],  c  is  the  sediment  concentration  [particle/(m.m.m)], 
and  r_s  is  the  mass  per  sediment  particle  [kg/particle].  Further: 


sources  —  sinks  =  .Dr(rf  l )  -f  D^(rfi)  —  &(T(rf  {)  —  |tjs(r,  /)|j  4-  /)  (5) 

T(r>  t)  =  9h( r,  t )  S( r)  -  rCT]  (6) 

where  T  is  the  sediment  transport  capacity  [kg/(ms)],  D_r  is  the  rill  erosion  or  deposition  rate,  D_i  is  the 
interrill  contribution  [kg/(m.m.s)],  C  is  the  first-order  reaction  coefficient  dependent  on  soil  and  cover 
[1/m],  r_w  is  the  mass  density  of  water  [kg/m.m.m],  S(r)= grad  z(r)  is  the  slope  [m/m],  z(r)  is  the 
elevation  [m],  Kt  is  the  transport  capacity  coefficient,  g=9.8l  is  gravitational  acceleration  [m/s.s.s],  We 
assume  that  the  critical  sheer  stress  tau_c  is  negligibly  small. 

Rill  detachment  and  deposition  are  proportional  to  the  difference  between  transport  capacity  and 
sediment  load  (eq.  5).  This  relationship  defining  the  interaction  between  sediment  load  and  transport 
capacity  (Foster  and  Meyer,  1972)  is  based  on  a  stream  power  concept  (Haan  1994)  and  can  be 
expressed  as: 

DJDc  +  q,/T=\  0) 

where  D_c-CT-K_d  r_w  gh.(.)S(.)  is  the  detachment  capacity  and  K_d  is  the  detachment  capacity 
coefficient  (rill  erodibility). 

We  solve  the  continuity  equation  in  2D  form  for  steady  state  water  flux  with  small  diffusive  term  with 
amplitude  d,  rewritten  in  the  operator  form  as: 

C^r)  =  ~|v2f(r)  +  V-  fe(r)v(r)5  +  Ce(r)iv(r)!  (8) 

and  then  the  erosion  equation  is 

Cq{  r)  =  CT(r)  +  D,-(r)  (9) 

where  rho=r_s  c(.)h(. ).  The  interpretation  of  the  equation  (8)  is  clear:  the  first  term  is  the  diffusion 
(which  in  our  case  is  very  small,  represents  the  smoothing  component  of  the  soil  transport),  the  second 
term  is  the  drift  driven  by  the  water  velocity  v(r)  and  the  third  term  is  the  ’potential’  which  is  dependent 
on  the  velocity  magnitude:  the  larger  the  velocity,  the  smaller  the  concentration  of  sediment. 

Net  erosion  and  deposition  is  then  estimated  as  a  divergence  of  sediment  flux.  Further  details  about  this 
approach  and  comparisons  with  previous  estimations  of  erosion/deposition  by  the  directional  derivative 
(Mitasova  et  al.  1995,  Wilson  and  Moore  1992)  will  be  given  elsewhere  Mitas  and  Mitasova,  (in  prep.). 

The  equation  (8)  is  solved  analogically  as  the  equation  for  water  flux  using  a  stochastic  method  (Mitas 
and  Mitasova,  in  prep.),  illustrated  by  Movie  l.: 
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Movie  1.  Solution  of  erosion  equations  by  Monte  Carlo,  illustrated  by  a  surface  representing  the 
sediment  flux  and  by  terrain  with  draped  erosion/deposition 


The  animation  shows  that  the  approximate  estimation  of  sediment  flux  is  reached  for  a  relatively  small 
number  of  Monte  Carlo  sampling  points.  Accurate  calculation  of  erosion/deposition  estimated  by 
derivatives  of  sediment  flux,  which  are  very  sensitive  to  statistical  noise,  can  be  estimated  to  a  given 
accuracy  either  by  employing  large  number  of  walkers  or  by  smoothing  out  the  noise  numerically. 


Influence  of  uniform  soil  and  cover  parameters 

In  the  erosion  model  the  influence  of  soil  and  cover  is  represented  by  the  following  basic  parameters: 
Mannings  n  ,  detachment  rate  coefficient  (erodibility)  Kd,  and  sediment  transport  coefficient  Kt.  These 
coefficients  are  functions  of  soil  and  cover  properties  such  as  soil  texture,  canopy,  roots,  management 
practices  etc.,  and  their  estimation  and  development  of  various  adjustment  factors  is  described  in  WEPP 
(Flanagan  and  Nearing  1995).  Constants  for  estimation  of  detachment  and  sediment  transport  capacity 
are  still  under  development  and  detailed  discussion  of  these  parameters  is  beyond  the  scope  of  this 
paper.  However,  we  will  use  the  following  examples  to  elucidate  the  role  of  these  parameters  in 
modeling  spatial  distribution  of  areas  with  erosion  or  deposition. 

The  most  important  parameter  controlling  the  borderline  between  the  erosion  and  deposition  areas  is  the 
first-order  reaction  coefficient  C  related  to  the  ratio  of  detachment  capacity  and  sediment  transport 
capacity.  In  the  first  example,  we  simulate  the  situation  when  the  study  area  has  constant  transport 
capacity  coefficient  Kt  but  detachment  capacity  coefficient  Kd  increases,  so  that  the  ratio  C  increases 
from  0.0005  to  100  (Movie  2). 
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Movie  2.  Change  in  the  spatial  distribution  of  sediment  flux  and  erosion/deposition  due  to  the  change  in 

C  with  increasing  Kd  and  Kt=  const. 


For  small  values  of  C  (e.g.,  clay  with  very  low  fall  velocities  and  low  detachment  capacity)  water  has 
the  power  to  carry  almost  all  detached  sediment  into  the  stream.  For  values  C>  1  ( e.g.,  sandy  soils  with 
relatively  high  fall  velocities  which  detach  easily),  the  sediment  flux  quickly  reaches  the  sediment 
transport  capacity  and  deposition  occurs  relatively  high  in  the  hillslope.  This  is  the  case  of  transport 
limited  erosion/deposition  modeled  by  a  simplified  approach  described  e.g.  by  Mitasova  et  al.  1995  and 
Moore  and  Wilson  1992.  For  both  cases  the  magnitude  of  sediment  flux  in  the  stream  remains  the  same 
while  the  distribution  of  erosion  and  deposition  over  the  landscape  changes  significantly.  This 
simulation  is  a  good  example  which  shows  that  calibrating  the  erosion  model  using  only  the  observed 
values  of  sediment  flux  at  an  outlet  does  not  guarantee  correct  predictions  of  erosion/deposition  on  the 
complex  hillslopes  within  the  watershed. 

Change  in  Kt  while  Kd  is  constant  also  changes  the  spatial  distribution  of  erosion/deposition,  however 
there  is  a  big  difference  in  the  amount  of  sediment  delivered  to  streams.  For  Kt«Kd  and  C»l  most  of 
the  detached  material  deposits  before  it  enters  the  stream,  for  Kt»Kd  and  C«l,  there  is  only  a  very 
small  deposition  and  most  of  the  detached  material  is  delivered  into  streams  (Movie  3).  This  example 
illustrates  how  the  potential  changes  in  soil  properties  and  cover  which  increase  transport  capacity  can 
trigger  severe  erosion. 


Movie  3.  Change  in  the  spatial  distribution  of  sediment  flux  and  erosion/deposition  due  to  the  change  in 

C  with  increasing  Kt  and  Kd-comt. 


If  C=const.  and  Kt,  Kd  change  with  the  same  rate,  the  spatial  distribution  of  erosion/deposition  is  the 
same,  and  only  its  magnitude  changes  (Figure  2.) 
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Figure  2.  Change  in  the  magnitude  of  erosion/deposition  with  C=1  and  increase  in  both  Kt  and  Kd; 


a)n=0. 1,  Kd=Kt==0.0003  (grass  on  sandy  soil);  b)n=0.05,  Kd=Kt=0.03  (bare  sandy  soil).  Surface 
topography  represents  the  sediment  flux,  color  is  the  erosion/deposition. 

It  is  important  to  note  that  the  parameters  Kt  and  Kd  which  are  dependent  on  the  soil  and  cover 
properties  are  interrelated  and  the  change  in  one  parameter  is  usually  accompanied  with  the  change  in 
the  second  parameter.  As  we  have  demonstrated,  it  is  their  ratio,  which  plays  important  role  in  the 
spatial  distribution  of  erosion  and  deposition.  With  better  understanding  of  the  physical  basis  of  these 
parameters  the  analysis  outlined  above  can  be  used  for  identification  of  those  soil  and  cover  properties 
which  can  be  targeted  for  the  most  effective  erosion  prevention. 

Erosion/deposition  for  spatio-temporal  changes  in  soil  and  cover  parameters 

Erosion  process  is  highly  dynamic  and  its  temporal  variability  can  be  modeled  at  various  time  scales 
from  minutes  (event  based  models  such  as  ANSWERS  or  AGNPS),  days  (WEPP),  to  geological  time 
(Moglen  and  Brass  1994).  For  land  use  management  applications  we  adapt  the  concept  used  in  WEPP 
and  we  simulate  erosion  under  steady  state  flow  for  variable  climatic,  soil  and  cover  parameters  as  they 
change  during  the  year. 

We  have  used  data  from  the  experimental  farm  of  the  Technical  University  in  Munchen  in  Germany 
(courtesy  Dr.  Karl  Auerswald)  such  as  elevation,  soil  core  data  and  current  land  use  and  we  simulated 
erosion/deposition  under  various  land  cover  and  climatic  conditions,  using  simulated  data  from  WEPP. 

Comparison  of  results  with  the  spatial  distribution  of  colluvial  deposits  (Figure  3)  and  with  the  pattern 
of  linear  erosion  after  a  150  year  storm  (Figures  3, 4)  indicates,  that  for  this  area  with  mostly  sandy  soils, 
the  terrain  controls  the  long  term  spatial  pattern  of  deposition  reflected  in  colluvial  deposits  observed  in 
mostly  concave  areas  (Figure  3a,  Movie  4).  Land  cover  has  more  significant  impact  on  the  magnitude  of 
erosion  and  short  term  linear  erosion  features  (Figure  4). 
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Figure  3.  Comparison  of  a)  observed  spatial  distribution  of  colluvial  deposits  (depth  in  cm)  and  linear 
erosion  features  after  the  150  year  storm  (red  lines)  with  b)  simulated  erosion/deposition  with 
homogeneous  bare  soil  conditions. 


Figure  4.  Simulated  b)  sediment  flux  and  b)  erosion/deposition  with  the  incorporation  of  the  a)  current 


land  use  influence. 


Movie  4.  Slicing  through  colluvial  deposits 

We  have  also  simulated  the  changes  of  erosion/deposition  during  the  year  due  to  the  changes  in  rainfall 
and  land  cover  illustrated  by  Movie  5.  and  Figure  5. 


Movie  5.  Change  in  the  land  cover  due  to  the  plant  growth  and  harvest 


Figure  5.  Impact  of  changes  in  rainfall  and  cover  on  erosion/deposition  distribution  and  sediment  flux  in 
May,  September  and  October. 

The  highest  risk  of  erosion  was  predicted  in  October,  when  there  is  minimum  cover  and  enough  rainfall 
to  produce  significant  runoff.  Although  May  had  the  most  intense  storm,  the  erosion  was  lower  due  to 
the  good  cover  provided  by  both  the  grass  area  and  the  agricultural  area  with  winter  wheat.  - 


LANDUSE  OPTIMIZATION 


Human  activity  changes  character  or  properties  of  landscape  components  (e.g.  cover)  which  can  be 
represented  as  an  appropriate  change  in  the  corresponding  field.  These  changes  influence  the  natural 
phenomena  through  various  interactions  which  can  be  described  by  the  governing  master  equations. 


Well-defined  and  quantified  impact  of  human  actions  on  nature  enable  researchers  to  formulate 
objectives  or  costs  in  order  to  either  predict  the  future  development  or,  more  importantly,  to  achieve  a 
desired  sustainable  development.  The  desired  objectives  can  include,  e.g.,  maximization  of  land  use  for 
production  or  military  training  with  minimized  impact  on  environment,  or  prevention  of  unacceptable 
changes  in  environment  in  the  given  time  horizon  with  minized  costs  (Johnston  and  Hopkins  1994). 
Because  of  the  extremely  complex  nature  of  the  problem,  the  optimization  tasks  are  often  out  of  the 
scope  of  ordinary  techniques  as  they  involve  multivariate  fields  (possibly  evolving  in  time)  and  also 
because  of  the  special  type  of  human  action  (like  instantaneous  point  sources  such  as  contamination 
which  spreads  out  in  the  time  horizon  of  a  few  years  or  clear  cut  of  a  forest  with  consequent  erosion). 
Therefore  this  problem  requires  a  formulation  of  general  methods  which  can  deal  with  complicated 
types  of  "configuration  or  state  spaces".  For  our  case  we  can  define  the  state  space  as  a  set  of  fields  (i.e., 
a  particular  set  of  multivariate  functions)  which  describe  components  of  the  studied  phenomenon. 
Available  information  and  models  such  as  initial  fields  values,  are  provided  by  a  GIS  and  are  used  as 
inputs  into  the  master  equations  and  their  solvers.  In  addition,  we  need  to  express  the  objective  (cost) 
functional  which  is  to  be  minimized  within  given  constraints.  The  constraints  can  be  formulated  in  the 
form  of  “external”  fixed  influences  (e.g.  part  of  land  cover  which  cannot  be  changed),  thresholds  on 
evolving  fields  (e.g.  erosion  beyond  certain  level  is  unacceptable)  and  so  forth.  The  general  form  of  the 
cost  functional  can  be  given  as: 

/  = J  j  drdtFdz^i)})  (10) 

where  /  is  the  cost  functional,  zj  are  the  input  spatio-temporal  fields,  F  is  a  function  which  determines 
of  cost  for  a  particular  set  of  z_i  (point  in  the  state  space).  In  general,  a  minimization  of  (10)  can  be  a 
very  complex  task.  In  order  to  carry  out  the  minimization  of  the  functional  (10)  we  have  to  define  the 
following: 

•  "distance"  in  the  state  space 

•  efficient  representation  of  the  fields  which  can  be  varied  by  the  human  impact  (e.g.  using 
appropriate  basis  function  expansions) 

•  "movement"  of  the  space  of  field  configurations 

Another  important  task  is  to  formulate  efficient  minimization  strategies.  Because  we  are  dealing  with  a 
multivariate  problem  which  often  involves  non-linearities  the  cost  functional  can  have  many  local 
minimas.  This  requires  use  of  robust  minimization  methods  such  as  simulated  annealing  or  genetic 
algorithms. 

A  simple  example  of  using  GIS  and  the  erosion  model  to  optimize  land  use  by  finding  a  more  effective 
spatial  distribution  of  protective  grass  cover  while  keeping  the  ratio  between  the  agricultural  area  and 
area  protected  by  grass  constant  is  illustrated  in  Figures  6  and  7.  Under  the  current  land  use  (Figure  6a), 
there  is  still  a  significant  amount  of  sediment  delivered  to  the  stream  (Figure  6b)  with  strong  potential  of 
creating  rills  and  gullies  (Figure  6c,  dark  red)  which  is  in  agreement  with  observations  of  big  storm 
-  effects,  presented  in  Figure  4.  Redesigning  the  land  use  so  that  the  protective  grass  cover  is  located  in 
the  highest  erosion  risk  areas  (Figure  7a)  can  dramatically  reduce  soil  loss  and  sediment  delivery  to  the 
streams  (Figure  7b,c).  The  crest  in  sediment  flux  in  areas  with  observed  gullies  disappears  and  is 
replaced  by  light  deposition  caused  by  the  decrease  in  water  velocity  in  the  grass  strip  (Figure  7c). 
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Figure  6.  Current  situation:  a)  land  use,  b)  simulated  sediment  flux,  c)  erosion/deposition 
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Figure  7.  Situation  with  the  grass  area  spatially  distributed  for  more  effective  erosion  protection  :  a) 
land  use,  b)  simulated  sediment  flux,  c)  erosion/deposition 

CONCLUSION 

We  have  presented  an  approach  to  modeling  of  landscape  processes  in  an  advanced  GIS  environment 
which  is  based  on  the  following  developments  : 

•  We  have  used  multivariate  regularized  splines  with  tension  for  scattered  data  interpolations  and 
continuous  field  representations  for  a  multitude  of  tasks  such  as  processing  of  input  data,  analysis 
and  presentation  of  simulation  results. 

•  We  have  developed  and  employed  Monte  Carlo  methods  for  solving  both  water  and  sediment  2D 
transport  problems.  This  approach  proved  to  be  robust  and  can  be  relatively  easily  generalized  to 
include  effects  currently  omitted  (e.g.,  3D  infiltration  process).  The  stochastic  methods  are  also 
very  well  suited  for  distributed  parallel  computing. 

•  Extending  the  erosion  model  from  traditional  ID  water  and  sediment  flow  equations  to  fully  2D 
fields,  supported  by  the  GIS  implementation,  provided  a  new  insight  into  the  functioning  of  these 
models  in  a  complex  realistic  landscape. 

•  Fully  integrated  visualization  based  on  multiple  dynamic  surfaces  helped  in  various  stages  of 
development,  evaluation  and  applications  of  complex  models,  where  interaction  between  the 
spatial  fields  is  important. 

•  On  the  basis  of  these  achievements  we  were  able  to  produce  landuse  scenario  with  a  potential  of 
significantly  decreased  erosion  in  a  fully  distributed  manner. 

Important  computational  components  of  our  approach  such  as  interpolation  tool,  water  flux  solver, 
sediment  flux  solver,  scenario  optimizer  (under  current  development)  are  built  as  functional  units  with 
well-defined  input,  output  and  controls.  Therefore  these  units  can  be  used  either  as  separate  tools  or 
within  open  GIS  frameworks  such  as  GRASS  depending  on  the  computational  environments  and  size  of 
tasks. 

Using  the  presented  application  as  an  example,  we  believe  that  the  GIS  in  future  can  become  not  only  a 


powerful  tool  for  providing  and  analyzing  spatial  information,  but  by  extending  its  capabilities  as  a 
simulation  and  optimization  tool,  it  can  allow  its  users  to  find  unexpected  solutions  of  land  management 
problems  leading  to  practices  which  can  be  more  effective  at  lower  cost  than  the  currently  known 
conservation  approaches. 
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